library(pander)
library(car)
ydat <- read.csv("yDat.csv", header=TRUE)
My final guess at your model.
\[ Y_i = \beta_0 + \beta_1 X_{10} + \beta_2 X_{10}^2 + \beta_3 X_{10}^3 + \beta_3 X_{10}^4 + \beta_3 X_{10}^5 + \beta_3 X_{10}^6 + \beta_3 X_{10}^7 + \beta_3 X_{10}^8 + \beta_3 X_{10}^9 + \beta_3 X_{10}^{10} + \beta_3 X_{10}^{11} + \beta_3 X_{10}^{12} + \epsilon_i \]
with my 95% confidence interval estimates of the coefficients as follows.
final.lm <- lm(Y ~ X10 + I(X10^2) + I(X10^3) + I(X10^4) + I(X10^5) + I(X10^6) + I(X10^7) + I(X10^8) + I(X10^9) + I(X10^10) + I(X10^11) + I(X10^12), data=ydat)
mytable <- round(confint(final.lm, level=1-0.05/length(final.lm$coef)), 2)
betas <- paste0("$\\beta_", 0:(length(final.lm$coef)-1), "$")
rownames(mytable) <- betas
colnames(mytable) <- c("Lower", "Upper")
pander(mytable)
| Lower | Upper | |
|---|---|---|
| \(\beta_0\) | 0 | 0 |
| \(\beta_1\) | 1 | 1 |
| \(\beta_2\) | 0 | 0 |
| \(\beta_3\) | 0.17 | 0.17 |
| \(\beta_4\) | 0 | 0 |
| \(\beta_5\) | 0.07 | 0.07 |
| \(\beta_6\) | 0 | 0 |
| \(\beta_7\) | 0.04 | 0.04 |
| \(\beta_8\) | 0 | 0 |
| \(\beta_9\) | 0.03 | 0.03 |
| \(\beta_10\) | 0 | 0 |
| \(\beta_11\) | 0.03 | 0.03 |
| \(\beta_12\) | 0.01 | 0.01 |
Glance at the data.
Certainly X10 matters.
pairs(ydat)
lm1 <- lm(Y ~ X10, data=ydat)
summary(lm1)
##
## Call:
## lm(formula = Y ~ X10, data = ydat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0054995 -0.0015366 -0.0001979 0.0012307 0.0135198
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0007937 0.0004862 1.632 0.111
## X10 1.0294709 0.0032153 320.177 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.003075 on 38 degrees of freedom
## Multiple R-squared: 0.9996, Adjusted R-squared: 0.9996
## F-statistic: 1.025e+05 on 1 and 38 DF, p-value: < 2.2e-16
plot(lm1, which=1)
Lots of multicollinearity here. I played with it for a while, and kept adding higher order polynomials until things broke with the final model I am presenting to you.
palette(colorRampPalette(c("gray","black"))(length(ydat$Y)))
pairs(cbind(R=lm1$res, ydat), col=as.factor(ydat$Y))
lm2 <- lm(Y ~ X10 + I(X10^2) + I(X10^3) + I(X10^4) + I(X10^5) + I(X10^6) + I(X10^7) + I(X10^8) + I(X10^9) + I(X10^10) + I(X10^11) + I(X10^12), data=ydat)
plot(Y ~ X10, data=ydat)
b <- coef(lm2)
curve(b[1] + b[2]*x + b[3]*x^2 + b[4]*x^3 + b[5]*x^4 + b[6]*x^7 + b[7]*x^8 + b[8]*x^9 + b[9]*x^10 + b[10]*x^11 + b[11]*x^12, add=TRUE)
plot(lm2, which=1)
makeTable9.3(ydat)
## Reordering variables and trying again:
## p model
## 1 1 (intercept)
## 2 2 X10
## 3 3 X10,X18
## 4 4 X7,X9,X10
## 5 5 X7,X9,X10,X18
## 6 6 X1,X7,X9,X10,X17
## 7 7 X1,X3,X7,X9,X10,X17
## 8 8 X1,X7,X9,X10,X17,X18,X20
## 9 9 X1,X3,X7,X9,X10,X17,X18,X20
## 10 10 X1,X3,X7,X9,X10,X12,X17,X18,X20
## 11 11 X1,X3,X4,X7,X9,X10,X12,X13,X17,X18
## 12 12 X1,X3,X4,X7,X9,X10,X12,X13,X15,X17,X18
## 13 13 X1,X4,X5,X7,X8,X10,X11,X12,X13,X16,X17,X18
## 14 14 X1,X4,X5,X7,X8,X10,X11,X12,X13,X16,X17,X18,X19
## 15 15 X1,X4,X5,X7,X8,X9,X10,X12,X15,X16,X17,X18,X19,X20
## 16 16 X1,X4,X6,X7,X8,X9,X10,X11,X12,X13,X14,X15,X17,X18,X19
## 17 17 X1,X3,X4,X5,X6,X7,X8,X9,X10,X12,X13,X14,X15,X17,X18,X19
## 18 18 X1,X3,X4,X6,X7,X8,X9,X10,X11,X12,X14,X15,X16,X17,X18,X19,X20
## 19 19 X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X12,X13,X14,X15,X16,X17,X18,X19
## SSEp R2p R2ap Cp AICp SBCp
## 1 9.697694e-01 0.0000000 0.0000000 NA NA NA
## 2 3.593445e-04 0.9996295 0.9996197 70.671131 -345.2893 -308.6435
## 3 2.527134e-04 0.9997394 0.9997253 41.017793 -357.3703 -319.0357
## 4 1.860036e-04 0.9998082 0.9997922 23.215026 -367.6299 -327.6064
## 5 1.639249e-04 0.9998310 0.9998116 18.660985 -370.6842 -328.9718
## 6 1.222478e-04 0.9998739 0.9998554 8.289154 -380.4185 -337.0173
## 7 1.137025e-04 0.9998828 0.9998614 7.752490 -381.3171 -336.2270
## 8 1.033109e-04 0.9998935 0.9998702 6.667761 -383.1508 -336.3718
## 9 9.473496e-05 0.9999023 0.9998771 6.122000 -384.6172 -336.1493
## 10 8.778793e-05 0.9999095 0.9998823 6.059782 -385.6636 -335.5068
## 11 8.304397e-05 0.9999144 0.9998848 6.651539 -385.8857 -334.0401
## 12 7.943735e-05 0.9999181 0.9998859 7.580916 -385.6618 -332.1272
## 13 7.507399e-05 0.9999226 0.9998882 8.285656 -385.9216 -330.6981
## 14 7.331834e-05 0.9999244 0.9998866 9.764493 -384.8681 -327.9558
## 15 7.119715e-05 0.9999266 0.9998855 11.134821 -384.0424 -325.4412
## 16 6.898303e-05 0.9999289 0.9998844 12.477561 -383.3061 -323.0160
## 17 6.524291e-05 0.9999327 0.9998859 13.367307 -383.5358 -321.5569
## 18 6.455524e-05 0.9999334 0.9998820 15.163175 -381.9597 -318.2919
## 19 6.400555e-05 0.9999340 0.9998774 17.000000 -380.3017 -314.9450
## PRESSp
## 1 NA
## 2 0.0006198257
## 3 0.0005619685
## 4 0.0006821237
## 5 0.0006165820
## 6 0.0007122717
## 7 0.0006675873
## 8 0.0006351166
## 9 0.0005732631
## 10 0.0005056526
## 11 0.0004775183
## 12 0.0004571693
## 13 0.0004289778
## 14 0.0004199173
## 15 0.0004162650
## 16 0.0005118786
## 17 0.0005226417
## 18 0.0005263321
## 19 0.0005660617